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We propose a new optimization method based on a demagnetization procedure well known in 
magnetism. We show how this procedure can be applied as a general tool to search for optimal 
solutions in any system where the configuration space is endowed with a suitable 'distance'. We test 
the new algorithm on frustrated magnetic models and the traveling salesman problem. We find that 
the new method successfully competes with similar basic algorithms such as simulated annealing. 
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It has been observed long time ago that disordered 
materials can be brought into a remarkably stable state 
through annealing, i.e. cooling down the material rather 
slowly. This simple observation inspired Kirkpatrick, 
Gelatt, and Vecchi in their pioneering work [0 to inves- 
tigate how close this annealing procedure takes models 
with glassy properties to their ground state, and lead 
them to the invention of the by now widely used sim- 
ulated annealing (SA) procedure. The SA revealed the 
crucial role the external noise can play in optimization: 
Thermal noise can help to escape high-energy local min- 
ima. In the present work we investigate another proce- 
dure that is commonly used to demagnetize disordered 
magnets and is experimentally known to result in a very 
stable state: The application of an oscillating external 
field (see Figure 1). This procedure makes use of an- 
other type of noise which is typical in magnetic systems, 
namely random external fields. As we shall demonstrate 
below for various models, a simple generalization of this 
zero-temperature 'a.c. demagnetization' is able to give 
systematically better and better approximations to the 
ground state of these models, and is in many cases 5-10- 
times faster than SA. We shall show how this method can 
be applied to practically any disordered model, thereby 
resulting in a new optimization procedure, that we call 
hysteretic optimization (HO). 

Finding optimal solutions of complex problems de- 
pending on a large number of parameters is an equally 
important and difficult task Examples range from 
integrated circuit design, through portfolio selection on 
the stock market [|| and calculating protein folding, to 
teaching artificial neural networks, to name a few. The 
simultaneous presence of randomness and frustration is 
what makes these problems so hard: Disorder is caused 
by the non-regular dependence of the quality of the so- 
lution on the configuration, and frustration is brought in 
by the competition of mutually exclusive different "good" 
properties. As a result, on one hand, a naiv search often 
gets stuck in spurious minima while, on the other hand, 
comparably good solutions can be found with quite dif- 



ferent configurations. 
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FIG. 1. The total internal energy (without the Zeeman 
term) of a three dimensional (L = 10) Edwards- Anderson 
spin glass during a.c. demagnetization. The inset shows the 
corresponding magnetization curve. 

It is important to note that most hysteretic systems 
fulfill the above requirements of complexity. Hysteresis 
implies the presence of many metastable states caused 
mainly by disorder. In the case of magnetic materials, 
the other important ingredient, frustration, is furnished 
by the magnetostatic interaction, which can be ferro- or 
antiferromagnetic depending on the relative orientation 
of the dipoles. This analogy between magnetic sytems 
and optimization suggests that part of knowledge ac- 
cumulated through decades by hysteresis research 
will eventually prove useful in optimization (and vice 
versa). Indeed, a simple but very frequently used hys- 
teresis model, the Preisach model , can be put in its 
ground state by a.c. demagnetization. 

In this letter, we proceed by studying in detail the en- 
ergetics of a.c. demagnetization for spin glass models. 
Then we show how this simple method can be modified 
to obtain a hysteretic optimization technique, which sys- 
tematically approaches the ground state , and we com- 
pare the new algorithm with SA. We also present how HO 
can be combined with cluster renormalization, and apply 
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the general formalism of HO for the traveling salesman 
problem (TSP). 

Many optimization problems can be formulated in 
terms of interacting Ising spins (a = ±1). Therefore 
we studied first two classical Ising spin glass models ||] 
with a Hamiltonian: 

JV N 

i,j i 

In the case of the Sherrington-Kirkpatrick (SK) model 
for any pair J tj = Zij/yN, where random 
gaussian number with zero mean and unit variance. In 
the other spin glass we considered, the three dimensional 
Edwards- Anderson (EA) model, Jij = Zij, but only for 
the nearest- neighbor couplings. The direction & = ±1 of 
the external field H randomly changes from site to site. 
Due to a spin-gauge symmetry, any choice of £j is equiv- 
alent, as long as it is not correlated with Jy. Using this 
symmetry, and also the analogy of the magnetic energy 
from equation ([[]), we call m = l/Nj^&o'i magnetiza- 
tion throughout the paper. These sytems are in the focus 
of research both in magnetism and in optimization. In 
magnetism their dynamics and the nature of the ordered 
state is still highly debated jj|,[L0|. From the optimiza- 
tion point of view, they are interesting because they are 
thought to be so-called NP-hard problems (ll]] . The ma- 
jor difference between the two models lies in the range 
of interactions. The short-ranged EA model may show 
some kind of clusterization ("droplets"), which is very 
unlikely for the infinite-ranged SK model. 
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FIG. 2. Size dependence of the integrated energy distri- 
bution of states produced by a.c. demagnetization for the SK 
model (7 = 0.9). E is the energy per spin measured from the 
ground state of each sample. 

We tested first the energetic properties of a.c. demag- 
netization on the SK model |n| . Many properties of this 
model can be calculated exactly J^] , and it is an ideal test 
system. We know, e.g., that its ground state energy is 
— 0.765/spin as N — > 00. Starting from a random initial 
state, a simple 'quench', when randomly chosen unsta- 
ble spins are aligned one-by-one with their local fields, 
hi = JijCj + H^i, ends up with (locally) stable states 
of energy around —0.70 per spin at H = 0. 



We started the a.c. demagnetization with an exter- 
nal field that saturates the spins, i.e. an H = Hi large 
enough that the Oi = & state (to = 1) be stable. Then 
we decreased the field H in small steps, each time mak- 
ing a few spins unstable, and 'quenching' the system to 
a nearby stable state using sequential dynamics. When 
H reached the 'turning point' H2 = — 7iiJi, we started 
to increase the field again until H3 = —72-^2, and turned 
back again and again at points H n = — 7 n _iif n _i, un- 
til the amplitude of the loop, \H n \, reached the order of 
the field steps. The resulting "spiraling" magnetization 
curve is shown in the inset of Figure [I]. We performed 
a.c. demagnetization with loops which were in average 
reduced by a factor of 7 = (7,;) = 0.9 after each turn- 
ing point (we used 7j-s between 0.8 and 1.0). We found 
that, as expected, lower 7 gives poorer results, but a 7 
even closer to one did not make a considerable improve- 
ment. However, randomly varying the turning points, 7$, 
increased the efficiency. 

Figure [2] shows for different system sizes the proba- 
bility of finding a configuration close to the ground state 
Jl3| of the SK model with a.c. demagnetization. For 
each N we averaged over several realizations of the cou- 
plings and ran a.c. demagnetization many times with 
different random field directions to generate the dis- 
tributions presented. The finite intercept at E = for 
smaller systems signals a finite probability of finding the 
ground state, but this probability goes to zero as the size 
increases: for very large systems a.c. demagnetization 
typically gives a state with an energy 0.012 x TV above 
the ground state. 
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FIG. 3. Fig. a: Probability density of finding a state of 
an SK model (TV = 1000) with energy E per spin above the 
ground state by a.c. demagnetization (symbols) and HO with 
1, 10, and 100 shake ups (solid lines). Fig.b: Average energy 
per spin measured from the ground state vs. running time for 
the SK model (N = 1000). 
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The probability density P(E) of finding a state with 
energy E produced by a.c. demagnetization for the 
N = fOOO SK and the f x 10 x 10 EA models is shown 
by the symbols in Figs. ||.a and [|. These measurements 
were performed for a single sample (instance). The a.c. 
demagnetization gives remarkable improvement with re- 
spect to the quench, but in both cases the final state is 
definitely above the ground state. We emphasize, how- 
ever, that a.c. demagnetization runs 3-10 times faster 
than SA (see Figs. ||.b and and can be a viable al- 
ternative in applications where a large number of opti- 
mizations is required, each of them on relatively smaller 
systems. 

As we mentioned earlier, unlike SA, increasing the fac- 
tor 7 closer to one, we only find a limited effect on the 
quality of solutions obtained by a.c. demagnetization. 
This is a serious problem, and finding the cause of it is 
a potentially important subject for future research. We 
speculate, however, that the answer is hidden in the ran- 
dom, but "frozen" correlations between our noise term 
£i and the generated low-lying states. These correlations 
may prevent us from reaching certain parts of the phase 
space. If so, all we have to do is to introduce different 
£,-s, since they have different correlations. But we have 
to do this without destroying the "good" correlations of 
previous £,-s. 

Maybe the simplest way to achieve this is "shaking up" 
the system: After finishing an a.c. demagnetization, H 
being zero, we are free to change £j to a different noise 
"direction" ^, while the final spin configuration of re- 
mains stable. Starting from this state, we increase the 
field in small steps to H m , quenching the spins after each 
step as before. From this point we start a new demagne- 
tization by decreasing the field to —"/iH m , and turning 
back e.t.c. At the and of the shake up we get a new 
state erf, and H is zero again, so we can choose a new 
and start another shake up. While subsequent shake ups 
with appropriate H m improve the solution on average, it 
is obviously helpful to get rid of "bad fluctuations", i.e. 
in case o\ had higher energy than of , we would start the 
new shake up from the latter again. It is this combina- 
tion of a.c. demagnetization and shake ups what we call 
hysteretic optimization (HO). 

The exact choice of the shaking amplitude H m is best 
obtained by testing it with a few discrete values: too 
small loops make hardly any change and too large field 
values take us back to a.c. demagnetization results. Our 
practice shows that the best choice for H m is around 
the coercive field, i.e. the field amplitude at which the 
magnetization reaches zero if started from saturation (~ 
0.3 on Fig. |). 

Figs. U and [| summarize our results of the shake-up 
HO calculations on the SK and EA models. Apparently, 
already the first few shake-ups give substantial improve- 
ment over a.c. demagnetization. Increasing further the 
number of shake-ups is less and less effective, but we 



get systematically closer and closer to the ground state, 
similar to SA. We would like to emphasize that apply- 
ing shake-ups after an a.c. demagnetization is straight- 
forward, requires very little extra coding and computer 
time. A systematic comparison of a.c. demagnetization, 
shake-up HO, and SA is shown in Figs. |].b and |[ 

In the case of the short-range EA model (inset of Fig- 
ure |]) shake-ups seem less effective, and for longer run- 
ning times SA reaches the efficiency of HO. However, typ- 
ical optimization problems are often closer to the long- 
range SK model where no particular dimensionality is 
present (e.g. stock market) ||. In this case (Fig. ^.b), 
HO appears to be better than SA even for long running 
times. 
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FIG. 4. Probability density of finding a state with energy 
E per spin above the ground state for the three-dimensional 
(L — 10) Edwards- Anderson model by a.c. demagnetization 
(symbols) and HO with 1, 10, and 100 shake ups (solid lines). 
The dashed line shows results of HO combined with RG. Inset: 
Typical energies as a function of running times for SA and HO 
with shake ups. 

We have to emphasize that, as well as SA the HO pre- 
sented in this paper should be considered as a building 
block, and it can be efficiently combined with other opti- 
mization methods, such as genetic algorith ms |14| , clus- 
ter rcnormalization group (RG) techniques [15|, or both 
fust . To demonstrate this, we combined HO with the 
RG in the following way: We generated clusters from 
n(~ 10 — 20) different low-energy configurations obtained 
by HO and using the clusterization technique of Ref. jl6| . 
Using these cluster spins simplifies the problem not just 
by decreasing the number of degrees of freedom, but also 
by removing dynamical barriers: turning over a bigger 
cluster through single spin flips might be energetically 
expensive because of the strong couplings inside the clus- 
ter. For the effective system of cluster spins HO can be 
applied again to generate n new (cluster) configurations 
and to define clusters of clusters, and so on. We used this 
relatively large n [jlTj on each RG level, and consequently 
many RG steps, to avoid the necessity of including a ge- 
netic algorithm. 

It is worth noting that there is a freedom in introducing 
the external field for the clusters, since the "magnetic 
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moment" /i of a cluster can be arbitrarily chosen. The 
simplest choice is /i = 1 for all clusters, but it might be 
useful to scale the local moments to insure that small and 
large clusters can compete with each other. The success 
of the RG approach for the EA model is due mostly to the 
short-range nature of the interactions (see Fig. We 
find that this method is able to reach the ground state 
in a time comparable to a recent very efficient genetic 
algorithm developed specifically for this problem 0,0 ■ 

Now we present a general scheme that makes it possi- 
ble to apply the present algorithm to practically any op- 
timization problem. In an optimization problem one has 
to minimize some cost function W(P), which is a map- 
ping from some configuration space {P} to real numbers. 
In order to apply HO to a system we need three impor- 
tant ingredients: (i) Dynamics. This is already sufficient 
to do SA, however, for the HO we need also a (ii) dis- 
tance d{P, Q) over the configuration space, and (iii) two 
reference states R±. These latter will be used to induce 
a random external field. 

We have a freedom in chosing these ingredients. How- 
ever, a 'good' dynamics must be such that an elementary 
step connects configurations with approximately equal 
cost functions, and its successive applications connect all 
configurations. Similarly, a 'good' distance correspond- 
ing to a given dynamics is such that distances do not 
change much during a single step. 

Having defined all these quantities we define the HO 
as the demagnetization of the Hamiltonian 

H(P) = W(P) + aHQ(aH)d(P, R a ), (2) 

a=± 

whith Q(x) being the step function. For the spin 
problems considered here the configurations are P — 
{<Jj}, the dynamics is single spin flip, the distance is 
d({ci}, = Yli \°i ~ an d the two reference states 
are R± = {±£l, .., ±fjv}- With these definitions Eq. (||) 
reduces to Eq. and the general HO becomes the one 
we applied to the spin problems. The basic role of the 
external field in Eq. (|^) is to force the system close to the 
states R± as H — > ±oo. 



Now we demonstarate how this method can be applied 
for the TSP. In this classical optimization problem the 
aim is to find the shortest path that goes through each 
of a given set of % = 1, .., N cities, visiting just once every 
one of them. In this problem the randomness is in the 
location Xj of the cities. The naive strategy of going 
always to the closest unvisited city fails, because typically 
there is a huge penalty for the last trip home — giving 
the major source of frustration. The configuration space 
is just given by all possible permutations P of the cities, 
and the cost function is simply the total length of the 
path, W(P) — 5Z i=1 |xp(i) — Xp(j +1 )|. As an elementary 
step we used the interchange of two cities jl9]], and a 
distance 

N 

d(P,Q) = 2|A i (P)-A i (Q)|, (3) 

i 

where Aj(P) = xp( j+1 j — Xp(j) denotes the vector con- 
necting the cities P(i) and P(i + 1). Our results for the 
traveling salesman are summarized in Fig. |^. We find 
a systematic improvement as the number of shake ups 
increases, and get better and better results in this case 
too. 

There is one more aspect of HO we would like to em- 
phasize. SA has the nice property of self-consistently 
telling us how the annealing schedule should be set: the 
specific heat is proportional to the fluctuations, so one 
might want to spend more time in the temperature re- 
gions of high specific heat. A similar measure of accuracy 
for HO might be the measured Preisach function of the 
system. This function can be obtained numerically from 
first order reversal curves ||, and once at hand, it can 
guide the annealing of the loops. We believe that this 
function has a lot of information about the metastable 
states of the system, although admittedly future research 
is needed to identify these features. 
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FIG. 5. Distribution of length generated by HO for a 
two-dimensional TSP with TV = 100. For comparison, we also 
show the best of 24 quenches, taking the same time as a.c. 
demagnetization . 
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